home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / POLYROOT.ZIP / DPRBM.FOR < prev    next >
Text File  |  1985-11-29  |  8KB  |  252 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE DPRBM
  5. C
  6. C        PURPOSE
  7. C           TO CALCULATE ALL REAL AND COMPLEX ROOTS OF A GIVEN
  8. C           POLYNOMIAL WITH REAL COEFFICIENTS.
  9. C
  10. C        USAGE
  11. C           CALL DPRBM (C,IC,RR,RC,POL,IR,IER)
  12. C
  13. C        DESCRIPTION OF PARAMETERS
  14. C           C      - DOUBLE PRECISION INPUT VECTOR CONTAINING THE
  15. C                    COEFFICIENTS OF THE GIVEN POLYNOMIAL. COEFFICIENTS
  16. C                    ARE ORDERED FROM LOW TO HIGH. ON RETURN COEFFI-
  17. C                    CIENTS ARE DIVIDED BY THE LAST NONZERO TERM.
  18. C           IC     - DIMENSION OF VECTORS C, RR, RC, AND POL.
  19. C           RR     - RESULTANT DOUBLE PRECISION VECTOR OF REAL PARTS
  20. C                    OF THE ROOTS.
  21. C           RC     - RESULTANT DOUBLE PRECISION VECTOR OF COMPLEX PARTS
  22. C                    OF THE ROOTS.
  23. C           POL    - RESULTANT DOUBLE PRECISION VECTOR OF COEFFICIENTS
  24. C                    OF THE POLYNOMIAL WITH CALCULATED ROOTS.
  25. C                    COEFFICIENTS ARE ORDERED FROM LOW TO HIGH (SEE
  26. C                    REMARK 4).
  27. C           IR     - OUTPUT VALUE SPECIFYING THE NUMBER OF CALCULATED
  28. C                    ROOTS. NORMALLY IR IS EQUAL TO IC-1.
  29. C           IER    - RESULTANT ERROR PARAMETER CODED AS FOLLOWS
  30. C                     IER=0  - NO ERROR,
  31. C                     IER=1  - SUBROUTINE DPQFB RECORDS POOR CONVERGENCE
  32. C                              AT SOME QUADRATIC FACTORIZATION WITHIN
  33. C                              100 ITERATION STEPS,
  34. C                     IER=2  - POLYNOMIAL IS DEGENERATE, I.E. ZERO OR
  35. C                              CONSTANT,
  36. C                              OR OVERFLOW IN NORMALIZATION OF GIVEN
  37. C                              POLYNOMIAL,
  38. C                     IER=3  - THE SUBROUTINE IS BYPASSED DUE TO
  39. C                              SUCCESSIVE ZERO DIVISORS OR OVERFLOWS
  40. C                              IN QUADRATIC FACTORIZATION OR DUE TO
  41. C                              COMPLETELY UNSATISFACTORY ACCURACY,
  42. C                     IER=-1 - CALCULATED COEFFICIENT VECTOR HAS LESS
  43. C                              THAN SIX CORRECT SIGNIFICANT DIGITS.
  44. C                              THIS REVEALS POOR ACCURACY OF CALCULATED
  45. C                              ROOTS.
  46. C
  47. C        REMARKS
  48. C           (1) REAL PARTS OF THE ROOTS ARE STORED IN RR(1) UP TO RR(IR)
  49. C               AND CORRESPONDING COMPLEX PARTS IN RC(1) UP TO RC(IR).
  50. C           (2) ERROR MESSAGE IER=1 INDICATES POOR CONVERGENCE WITHIN
  51. C               100 ITERATION STEPS AT SOME QUADRATIC FACTORIZATION
  52. C               PERFORMED BY SUBROUTINE DPQFB.
  53. C           (3) NO ACTION BESIDES ERROR MESSAGE IER=2 IN CASE OF A ZERO
  54. C               OR CONSTANT POLYNOMIAL. THE SAME ERROR MESSAGE IS GIVEN
  55. C               IN CASE OF AN OVERFLOW IN NORMALIZATION OF GIVEN
  56. C               POLYNOMIAL.
  57. C           (4) ERROR MESSAGE IER=3 INDICATES SUCCESSIVE ZERO DIVISORS
  58. C               OR OVERFLOWS OR COMPLETELY UNSATISFACTORY ACCURACY AT
  59. C               ANY QUADRATIC FACTORIZATION PERFORMED BY
  60. C               SUBROUTINE DPQFB. IN THIS CASE CALCULATION IS BYPASSED.
  61. C               IR RECORDS THE NUMBER OF CALCULATED ROOTS.
  62. C               POL(1),...,POL(J-IR) ARE THE COEFFICIENTS OF THE
  63. C               REMAINING POLYNOMIAL, WHERE J IS THE ACTUAL NUMBER OF
  64. C               COEFFICIENTS IN VECTOR C (NORMALLY J=IC).
  65. C           (5) IF CALCULATED COEFFICIENT VECTOR HAS LESS THAN SIX
  66. C               CORRECT SIGNIFICANT DIGITS THOUGH ALL QUADRATIC
  67. C               FACTORIZATIONS SHOWED SATISFACTORY ACCURACY, THE ERROR
  68. C               MESSAGE IER=-1 IS GIVEN.
  69. C           (6) THE FINAL COMPARISON BETWEEN GIVEN AND CALCULATED
  70. C               COEFFICIENT VECTOR IS PERFORMED ONLY IF ALL ROOTS HAVE
  71. C               BEEN CALCULATED. IN THIS CASE THE NUMBER OF ROOTS IR IS
  72. C               EQUAL TO THE ACTUAL DEGREE OF THE POLYNOMIAL (NORMALLY
  73. C               IR=IC-1). THE MAXIMAL RELATIVE ERROR OF THE COEFFICIENT
  74. C               VECTOR IS RECORDED IN RR(IR+1).
  75. C
  76. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  77. C           SUBROUTINE DPQFB    QUADRATIC FACTORIZATION OF A POLYNOMIAL
  78. C                               BY BAIRSTOW ITERATION.
  79. C
  80. C        METHOD
  81. C           THE ROOTS OF THE POLYNOMIAL ARE CALCULATED BY MEANS OF
  82. C           SUCCESSIVE QUADRATIC FACTORIZATION PERFORMED BY BAIRSTOW
  83. C           ITERATION. X**2 IS USED AS INITIAL GUESS FOR THE FIRST
  84. C           QUADRATIC FACTOR, AND FURTHER EACH CALCULATED QUADRATIC
  85. C           FACTOR IS USED AS INITIAL GUESS FOR THE NEXT ONE. AFTER
  86. C           COMPUTATION OF ALL ROOTS THE COEFFICIENT VECTOR IS
  87. C           CALCULATED AND COMPARED WITH THE GIVEN ONE.
  88. C           FOR REFERENCE, SEE J. H. WILKINSON, THE EVALUATION OF THE
  89. C           ZEROS OF ILL-CONDITIONED POLYNOMIALS (PART ONE AND TWO),
  90. C           NUMERISCHE MATHEMATIK, VOL.1 (1959), PP.150-180.
  91. C
  92. C     ..................................................................
  93. C
  94.       SUBROUTINE DPRBM(C,IC,RR,RC,POL,IR,IER)
  95. C
  96. C
  97.       DIMENSION C(1),RR(1),RC(1),POL(1),Q(4)
  98.       DOUBLE PRECISION C,RR,RC,POL,Q,EPS,A,B,H,Q1,Q2
  99. C
  100. C        TEST ON LEADING ZERO COEFFICIENTS
  101.       EPS=1.D-6
  102.       LIM=100
  103.       IR=IC+1
  104.     1 IR=IR-1
  105.       IF(IR-1)42,42,2
  106.     2 IF(C(IR))3,1,3
  107. C
  108. C        WORK UP ZERO ROOTS AND NORMALIZE REMAINING POLYNOMIAL
  109.     3 IER=0
  110.       J=IR
  111.       L=0
  112.       A=C(IR)
  113.       DO 8 I=1,IR
  114.       IF(L)4,4,7
  115.     4 IF(C(I))6,5,6
  116.     5 RR(I)=0.D0
  117.       RC(I)=0.D0
  118.       POL(J)=0.D0
  119.       J=J-1
  120.       GO TO 8
  121.     6 L=1
  122.       IST=I
  123.       J=0
  124.     7 J=J+1
  125.       C(I)=C(I)/A
  126.       POL(J)=C(I)
  127.       CALL OVERFL(N)
  128.       IF(N-2)42,8,8
  129.     8 CONTINUE
  130. C
  131. C        START BAIRSTOW ITERATION
  132.       Q1=0.D0
  133.       Q2=0.D0
  134.     9 IF(J-2)33,10,14
  135. C
  136. C        DEGREE OF RESTPOLYNOMIAL IS EQUAL TO ONE
  137.    10 A=POL(1)
  138.       RR(IST)=-A
  139.       RC(IST)=0.D0
  140.       IR=IR-1
  141.       Q2=0.D0
  142.       IF(IR-1)13,13,11
  143.    11 DO 12 I=2,IR
  144.       Q1=Q2
  145.       Q2=POL(I+1)
  146.    12 POL(I)=A*Q2+Q1
  147.    13 POL(IR+1)=A+Q2
  148.       GO TO 34
  149. C        THIS IS BRANCH TO COMPARISON OF COEFFICIENT VECTORS C AND POL
  150. C
  151. C        DEGREE OF RESTPOLYNOMIAL IS GREATER THAN ONE
  152.    14 DO 22 L=1,10
  153.       N=1
  154.    15 Q(1)=Q1
  155.       Q(2)=Q2
  156.       CALL DPQFB(POL,J,Q,LIM,I)
  157.       IF(I)16,24,23
  158.    16 IF(Q1)18,17,18
  159.    17 IF(Q2)18,21,18
  160.    18 GO TO (19,20,19,21),N
  161.    19 Q1=-Q1
  162.       N=N+1
  163.       GO TO 15
  164.    20 Q2=-Q2
  165.       N=N+1
  166.       GO TO 15
  167.    21 Q1=1.D0+Q1
  168.    22 Q2=1.D0-Q2
  169. C
  170. C        ERROR EXIT DUE TO UNSATISFACTORY RESULTS OF FACTORIZATION
  171.       IER=3
  172.       IR=IR-J
  173.       RETURN
  174. C
  175. C        WORK UP RESULTS OF QUADRATIC FACTORIZATION
  176.    23 IER=1
  177.    24 Q1=Q(1)
  178.       Q2=Q(2)
  179. C
  180. C        PERFORM DIVISION OF FACTORIZED POLYNOMIAL BY QUADRATIC FACTOR
  181.       B=0.D0
  182.       A=0.D0
  183.       I=J
  184.    25 H=-Q1*B-Q2*A+POL(I)
  185.       POL(I)=B
  186.       B=A
  187.       A=H
  188.       I=I-1
  189.       IF(I-2)26,26,25
  190.    26 POL(2)=B
  191.       POL(1)=A
  192. C
  193. C        MULTIPLY POLYNOMIAL WITH CALCULATED ROOTS BY QUADRATIC FACTOR
  194.       L=IR-1
  195.       IF(J-L)27,27,29
  196.    27 DO 28 I=J,L
  197.    28 POL(I-1)=POL(I-1)+POL(I)*Q2+POL(I+1)*Q1
  198.    29 POL(L)=POL(L)+POL(L+1)*Q2+Q1
  199.       POL(IR)=POL(IR)+Q2
  200. C
  201. C        CALCULATE ROOT-PAIR FROM QUADRATIC FACTOR X*X+Q2*X+Q1
  202.       H=-.5D0*Q2
  203.       A=H*H-Q1
  204.       B=DSQRT(DABS(A))
  205.       IF(A)30,30,31
  206.    30 RR(IST)=H
  207.       RC(IST)=B
  208.       IST=IST+1
  209.       RR(IST)=H
  210.       RC(IST)=-B
  211.       GO TO 32
  212.    31 B=H+DSIGN(B,H)
  213.       RR(IST)=Q1/B
  214.       RC(IST)=0.D0
  215.       IST=IST+1
  216.       RR(IST)=B
  217.       RC(IST)=0.D0
  218.    32 IST=IST+1
  219.       J=J-2
  220.       GO TO 9
  221. C
  222. C        SHIFT BACK ELEMENTS OF POL BY 1 AND COMPARE VECTORS POL AND C
  223.    33 IR=IR-1
  224.    34 A=0.D0
  225.       DO 38 I=1,IR
  226.       Q1=C(I)
  227.       Q2=POL(I+1)
  228.       POL(I)=Q2
  229.       IF(Q1)35,36,35
  230.    35 Q2=(Q1-Q2)/Q1
  231.    36 Q2=DABS(Q2)
  232.       IF(Q2-A)38,38,37
  233.    37 A=Q2
  234.    38 CONTINUE
  235.       I=IR+1
  236.       POL(I)=1.D0
  237.       RR(I)=A
  238.       RC(I)=0.D0
  239.       IF(IER)39,39,41
  240.    39 IF(A-EPS)41,41,40
  241. C
  242. C        WARNING DUE TO POOR ACCURACY OF CALCULATED COEFFICIENT VECTOR
  243.    40 IER=-1
  244.    41 RETURN
  245. C
  246. C        ERROR EXIT DUE TO DEGENERATE POLYNOMIAL OR OVERFLOW IN
  247. C        NORMALIZATION
  248.    42 IER=2
  249.       IR=0
  250.       RETURN
  251.       END
  252.